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We study a two-dimensional system of spin-polarized fermions on the kagome lattice at filling 
fraction / = 1/3 interacting through a nearest-neighbor interaction V . Above a critical interaction 
strength Vc a charge-density wave with a broken Z^. symmetry is stabilized. Using the unrestricted 
mean-field approximation, we present several arguments showing that elementary topological point 
defects in the order parameter bind a fractional charge. Our analysis makes use of two appealing 
properties of the model: (i) For weak interaction, the low-energy degrees of freedom are described 
by Dirac fermions coupled to a complex- valued mass field (order parameter), (ii) The nearest- 
neighbor interaction is geometrically frustrated at filling / = 1/3. Both properties offer a route to 
fractionalization and yield a consistent value ±1/2 for the fractional charge as long as the symmetry 
between the up and the down triangles of the kagome lattice is preserved. If this symmetry is 
violated, the value of the bound charge varies continuously with the strength of the symmetry- 
breaking term in the model. In addition, we have numerically computed the confining potential 
between two fractionally charged defects. We find that it grows linearly at large distances but can 
show a minimum at a finite separation for intermediate interactions. This indicates that the polaron 
state, formed upon doping the charge-density wave, can be viewed as a bound state of two defects. 



I. INTRODUCTION 

The concept of fermion fractionalization has been 
applied to a variety of condensed matter systems. 
Prominent examples are spin-charge separation in 
polyacetylene,^ fractionally charged excitations in the 
fractional quantum Hall states^ and magnetic monopoles 
in spin-ice.'^ In ah these examples, excitations carrying 
fractional quantum numbers with respect to the elemen- 
tary particles forming the system were found. The term 
fractionalization is stringently used only if well-defined 
excitations with fractional quantum number exist on all 
length scales. In gapped insulating systems in dimen- 
sions d > 2 this requires topological order of the ground 
state, ^ i.e. a ground state degeneracy which depends on 
the topology of the underlying system. Naturally, the 
effective low-energy theory describing the fractional ex- 
citations is a gauge theory in the deconfining phase, such 
as the Chern-Simons theory for the two dimensional frac- 
tional quantum Hall state or the Coulomb gauge theory 
for the three dimensional spin-ice materials and other 
frustrated magnetic systems. Electron fractionaliza- 
tion has also been discussed in the context of the high-T^ 
cuprates on the basis of the Z2 gauge theory.* In this 
article, however, we would like to use the concept of frac- 
tionalization in a less stringent way. Particularly, we are 
interested in phenomena where fractionalization occurs 
only up to a certain length scale which usually depends 
on temperatures and model parameters. This more gen- 
eral point of view allows one to cover a broader range 
of phenomena and to access this fascinating phenomena 
from distinct theoretical view-points. 

Several authors have stressed the field-theoretical point 
of view where fractionalized quantum numbers are car- 
ried by solitonic solutions of field theories which support 
isolated mid-gap states. ^'^^ In condensed matter physics, 
this route to fractionalization is well appreciated for the 



one-dimensional example discussed by Su, Schrieffer and 
Heeger where the soliton describes a domain wall sep- 
arating two degenerate dimerized ground states.^ More 
recently, a body of work has appeared^^"^° on general- 
izations of this concept to two dimensions and it has 
been argued that topological defects in the "kekule" or- 
der parameter offers an example for the solitonic frac- 
tionalization in two-dimensional graphene.^^ However, in 
contrast to the one-dimensional case, the energy cost as- 
sociated with a vortex in the complex-valued Bose field 
is not finite but grows with system size. In the contin- 
uum limit, the interaction between two vortices depends 
logarithmically on the distance and vortices can prolifer- 
ate above the Kosterlitz-Thouless temperature. On the 
lattice, however, they are always confined at sufficiently 
long distances. It was also pointed out that in order 
to heal the vortex at long distances, a coupling to an axial 
gauge field (which itself supports a vortex) can be intro- 
duced, thereby rendering the vortex energy finite. ^^'^^ 

Another perspective on the fractionalization phenom- 
ena emerges from models describing strongly interacting 
particles on a lattice with geometrical frustration. ^^'^^ 
In this class of models, the strong interaction enforces a 
local constraint and it is the violation of this local con- 
straint which carries a fractional charge. Clearly, there 
is a close relation to frustrated spin models and the afore 
mentioned spin-ice system is a prominent representative. 
In many cases, the frustrated particle interaction or spin 
exchange can be mapped on an effective hardcore dimer 
model. Removing one dimer introduces two monomers 
which, under certain circumstances, are well defined frac- 
tionalized excitations.^''"^® 

In this article, we want to make contact to both routes 
to fractionalization by studying topological point defects 
in a model of spin-polarized fermions on the kagome lat- 
tice subject to a nearest-neighbor interaction V. Pre- 
vious investigations^^'^* of this model at filling fraction 
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/ = 1/3 suggest a zero-temperature phase transition at 
a critical interaction Vc between the semi-metalhc Dirac 
hquid for V < Vc and a gapped and charge ordered 
state with a y/3 x y/3 reconstruction of the unit cell for 
V > Vc- Our discussion of topological point defects in 
the order parameter of the charge-density wave will make 
use of two important properties of the model. First, the 
low-energy degrees of freedom in the weakly interacting 
limit are well-described by Dirac-fermions coupled to the 
complex-valued order parameter which enters as a mass 
field. This offers the possibility for the solitonic fraction- 
alization mechanism in two-dimensions in analogy with 
graphene^^"^**'^^ and related systems. Second, the 
nearest-neighbor interaction is geometrically frustrated 
and the classical charge configurations with lowest en- 
ergy satisfy the "triangle rule". This constraint states 
that there is exactly one particle on every triangle of 
the kagome lattice. A local violation binds a fractional 
charge. Indeed, this possibility has recently been ex- 
plored in the strongly interacting limit using exact di- 
agonalization techniques and it has been argued that the 
defects carry a fractional charge ±1/2 and are asymptot- 
ically free in the large V limit. 

Our paper is organized as follows. In Sec. II we intro- 
duce the model and discuss the triangle rule and the effect 
of its local violation. In Sec. Ill we introduce the un- 
restricted mean-field approximation, discuss the leading 
instability of the Dirac liquid towards the charge-density 
wave and introduce a Ginzburg-Landau expansion of the 
free energy. This sets the stage for introducing topologi- 
cal point defects in Sec. IV and in Sec. V we numerically 
study solutions with point defects, compute the value 
of the bound charge and the confinement potential be- 
tween two defects. Eventually, in Sec. VI, we consider the 
weakly ordered state and establish a description in terms 
of Dirac fermions coupled to a complex- valued mass term. 



II. MODEL FOR CHARGE- ORDERED 
KAGOME LATTICE 

Our starting point is a tight-binding model of spin- 
polarized fermions on the kagome lattice at filling fraction 
/ = 1/3 subject to a nearest-neighbor repulsion V. The 
Hamiltonian is given by 
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Here, c^P annihilates (creates) a spin-polarized fermion 
on site i and rii = c^^c^. The hopping integral is denoted 
by t > and V > Q specifies the nearest-neighbor in- 
teraction. At several places in this article, we will also 
consider a term which enhances the hopping on the up 
triangles 
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FIG. 1. (a) The unit cell of the kagome lattice contains three 
sites labeled by 1,2 and 3. The vector r is an element of 
the underlying triangle lattice and denotes the location of the 
center of the up-triangles. The unit cell vectors are denoted 
by Oi and 02. (b) Density distribution in the charge density 
wave phase found in the model defined by Eq. (1). There is 
a -\/3 X -\/3 reconstruction of the unit cell. 



This term induces a bond-order wave which breaks the 
symmetry between the up and the down triangles. How- 
ever, unless otherwise stated, we set 5t = 0. The non- 
interacting {V — 0) band structure is obtained by diago- 
nalizing the matrix Ho{K) = —2tr{K) where 



r(K) 



cos(is:i/2) cos(i^2/2y 
cos(is:i/2) cos(K3/2) 
cos(fs:2/2) cos(if3/2) 



(3) 



Above, we have introduced K^, ^ K ■ a^j and (y — 
1, 2, 3) are given in Tab. I. There is a flat band at energy 
2t as well as two dispersing bands. It is well-known that 
at / = 1/3, the linearized band structure near the Fermi 
energy is described by two Dirac cones, similar to the 
situation found in graphene. 



A. Triangle rule in the atomic limit 

Let us now look at the atomic limit t = Q. It is known 
that at filling fraction 1/3 the interaction energy can be 
minimized by a macroscopic number of classical charge 
configurations. This fact becomes clear when rewriting 
the interaction Hamiltonian in real space as a sum over 
all triangles 5 of the kagome lattice in the following way: 

~ 2" 
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where Ng — J2ies denotes the total charge operator on 
the triangle 6 {d can label both an up or a down triangle) . 
Clearly, the interaction is lowest for configurations which 
fulfill the local constraint Ns — 1. Taking r to be the 
center of an up-triangle and using the labeling convention 
introduced in Fig. 1(a) and Tab. I this constraint takes 
the form 

ni{r) + n2{r) + nsi^r) ^ 1, (5a) 
nsi^r) + ni{r + a2) + n2ir + aa) = 1. (5b) 
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The first equation is written for the up and the second 
one for the down triangles. Thus, Hy is minimized by 
classical charge configurations with exactly one fermion 
on every triangle and the ground state is macroscopically 
degenerate. In the following, we refer to the local con- 
straint Eq. (5) as the "triangle rule"."^*^ It is instructive 
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TABLE I. Definitions of the lattice vectors used in this paper 
(in units of the lattice constant o) and the values of the in- 
ner product with the uniform ordering vector G = (87r/3,0) 
modulo 2tt. 




FIG. 2. (Color online.) A classical charge configuration which 
locally violates the triangle rule on the shaded (yellow) down 
triangle. We introduce a local charge density defined on every 
triangle as half the value of the sum of the charges on that 
triangle. In this way we see that the violation of the triangle 
rule carries a fractional charge 1/2. 



to see how the macroscopic degeneracy shows up in re- 
ciprocal space. Introducing the Fourier components 



we can write the interaction as 
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Here, N denotes the number of unit cells and we have 
introduced the vector notation 



(Q)t ^ n,{-Q) n^i-Q) n^i-Q] 
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The matrix r(Q) is given in Eq. (3). Its lowest eigen- 
value is equal to —1, independent of Q. It follows that 
the interaction energy is minimized by all charge config- 
urations which have Fourier components lying in the flat 
band. Indeed, it is straightforward to show that charge 
configurations which are proportional to the eigenvectors 
of the flat band fulfill the two constraints Eq. (5) in mo- 
mentum space. For Q = the conditions (5a) and (5b) 
are equivalent which is a manifestation of the quadratic 
band touching point at Q = in T{Q). 

The macroscopic degeneracy of the classical charge 
configurations is lifted for finite t. In particular, for 
t/V <^ 1, the model Eq. (1) can be mapped onto a quan- 
tum dimer model by identifying an occupied site of the 
kagome lattice with a dimer on the hexagonal lattice. ^^'^^ 
Thereby, ring exchange processes of order /V"^ in the 
original model translate into dimer flips in the dimer 
model which stabilizes a valence bond crystal with a 
\/3 X >/3 reconstructed unit cell.'^^ The kinetic energy 
gained by resonating plaquettes favors charge configu- 
rations which are connected by local dimer flips. The 
classical configuration which has most flippable plaque- 
ttes corresponds to the mean-field charge-density wave 
shown in Fig. 1(b). Note also that the constraint (5) 
maps onto a hardcore constraint for dimer coverings. 



B. Violation of triangle rule and fractional charge 

Let us now consider a classical charge configuration 
which locally violates the triangle rule Eq. (5) either for 
an up or a down triangle. An example is shown in Fig. 2 
where the triangle rule is violated on a single down trian- 
gle. In such a situation, the total charge density per unit 
cell depends on how it is measured! For example, if we 
measured it by summing the charges on the up triangles, 
we would conclude that there is exactly one particle in 
every unit cell. On the other hand, if we measured it 
by summing the charges on the down triangles, we would 
conclude that there is one particle missing in the unit cell 
which contains the empty down triangle. A more sensi- 
tive way which avoids this ambiguity is to introduce a 
charge density defined on every triangle as half the value 
of the sum of the charges on that triangle. This charge 
density is then defined on the hexagonal lattice formed 
by the center points of the triangles, see Fig. 2. If the 
triangle rule is fulfilled everywhere, there is a charge den- 
sity — 1/2 on every site of the hexagonal lattice (we as- 
sociate a charge —1 with a single particle). In this way 
we see that an empty triangle carries a fractional charge 
1/2 compared to a configuration which satisfies the tri- 
angle rule. Likewise, a triangle with two particles carries 
a charge —1/2 and with three particles a charge —1. 



III. MEAN-FIELD DESCRIPTION 

In this section we start with a conventional mean-field 
theory and discuss some properties of the uniformly or- 
dered system. In the Hartree approximation the density- 
density interaction is decoupled in the following way: 

n^{r)p^{r') + n^{r')p^{r) - p^{r)p^{r'). (9) 
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Here, we have introduced the expectation values 

p.(r) = K(r)) = l^p,(Q)e*Q-('^+'-^) (10) 
Q 

of the local densities with Fourier components Pv{Q)- 
In addition to the Hartree terms in Eq. (9) also the 
Fock terms have been considered in Ref. 27 for uni- 
form solutions. These terms tend to stabilize the semi- 
metallic phase and the critical interaction strength for 
the phase transition is V^^ ~ 3t, in good agreement 
with other methods. On the other hand, when keeping 
only the Hartree terms as in Eq. (9), the critical inter- 
action strength is smaller, w 2.2t. However, except 
for this shift, the qualitative aspects of the Hartree so- 
lution seems to be the same and to keep it simple, we 
use the decoupling Eq. (9). We note here that the sit- 
uation for filling fraction / = 2/3 is quite different be- 
cause complex Fock terms stabilize an interaction-driven 
topological insulator for arbitrary small nearest-neighbor 
interactions,^^ similar to what is found on the decorated 
honeycomb lattice at half filling. '^^ 

A. Mean-field triangle rule 

Any mean-field state is characterized by a self- 
consistent charge distribution {p^ir)} and the configura- 
tions with lowest energies fulfill the triangle-rule Eq. (5) 
on average. In Fourier space, we can write it for Q 7^ 
as 

= pi(Q)e-*«^/2 ^ p^(^Qy~^Q.n + p3(Q), (iia) 
= Pi(Q)e*^^/' + P2(g)e'^^/' + P3(Q), (lib) 

where Q^, — Q ■ a^, as before. If we introduce the vec- 
tor p{Q) — [pi{Q) P2{Q) PsiQ)]^ above condition is 
equivalent to the "flat-band" condition 

T(Q)p{Q) = -AQ)- (12) 
The mean-field interaction can then be written as 

Q,u Q,v 

where p{Q) satisfies Eq. (12). 

B. Leading instability 

In order to find the leading instability of the interact- 
ing system as function of the interaction V ^ we can study 
the static mean-field susceptibility associated with "flat- 
band" configurations [configurations which are compati- 
ble with Eq. (12)]: 

Xmf(Q) = r^l-ir- (13) 
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FIG. 3. (Color online.) The static susceptibility xq of the 
noninteracting system at T = associated with charge config- 
urations satisfying the triangle rule on average. Xq is largest 
at the corners of the hexagon forming the first Brillouin zone. 
The leading instability is a charge-density wave with one of 
the three ordering vectors Gi, G2 or G3. 

The leading instability at a fixed temperatures occurs 
at an ordering vector G which satisfies Vx% — 1 for 
the smallest value of V . This also defines the critical 
interaction = 1/xg- In Eq. (13) we have introduced 
the response function of the non-interacting system which 
probes charge configurations satisfying the triangle rule 
on average: 

Xq(T) = E [Go(zc.„)A(Q)Go(*c.„)A(-Q)] . 

iuJn 

Here, Go(ia;„) = l/[jwTi — + p\ is the Matsubara 
Green's function operator of the non-interacting system 
and LOn are fermionic Matsubara frequencies. The oper- 
ator A{Q) = e(g) • n{Q) with V{Q)e{Q) = -e{Q) and 
e{Q)-e{Q) = 1 enforces the triangle rule on average. The 
trace involves summation over K and the three bands of 
the noninteracting system. Figure 3 shows Xq{T = 0). 
The static susceptibility is largest at the corners of the 
hexagonal Brillouin zone and for the critical interaction 
at T = we find the numerical value V°{T = 0) w 2.33i. 
The leading instability therefore occurs at one of the 
three ordering vectors Gi, G2 or G3 connecting oppo- 
site corners of the hexagon. 

As mentioned earlier, this instability is a charge- 
density wave with a \/3 x \/3 reconstruction of the unit 
cell and is shown in Fig. 1(b). Instead of working with 
three different ordering vectors, we fix G = Gi and al- 
low for a complex phase of the charge-density wave order 
parameter, see below. In the following we choose 

G = K+ K (14) 

where K± denote the location of the two inequivalent 
Dirac points in the first Brillouin zone, 

K,=±(g,0). (15) 

Obviously, G couples the two Dirac points and it is this 
"nesting" which opens a gap above a critical interaction. 
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C. Ginzburg-Landau expansion 

From the triangle rule (11) it follows that pi{G) = 
P2{G) — /03(G) and we define the complex- valued order 
parameter of the charge-density wave as 

A = |A|e'^ = -^[pi(G) + P2(G) + p3(G)]. (16) 

Note that there is an overall phase freedom in the defi- 
nition of A and that the amplitude satisfies |A| < 2V/3. 
With the definition Eq. (16) of the order parameter, the 
interaction Hamiltonian for the uniform charge-density 
wave reduces to 



A' 



E 



3N 
2V^ 



A| 



(17) 



The free energy of a slowly varying charge-density wave 
in the continuum limit is given by the perturbative 
Ginzburg-Landau expansion in terms of A and its gradi- 
ent VA: 



F< 



CDW 



Fo 



dx dy 



a(F,T)|A|2-t-r;|A*VA| + «;|VA| 



-h7|Apcos(3v3) + ^|A| 



(18) 



Here, A = \^a'^/2 is the unit cell area of the kagome 
lattice and the coefficient /3 > stabilizes this expansion 
to order |A|^. The coefficient a{V,T) changes sign as 
function of the interaction strength or temperature and 
is given by 



a{V,T) 



1 



1 



V V^{T) 



V^{T) = 1/xUT). (19) 



The term proportional to cos(3(^) in the expansion 
Eq. (18) introduces an anisotropy as a result of the three- 
fold rotation symmetry of the triangular Bravais lat- 
tice. The numerical value of its prefactor at T = is 
7 This term acts as a pinning potential for 

the complex phase tp of the order parameter A and in 
the ground state, it assumes one of the three values 



-7r/3, ipB = tt/S, (fie 



(20) 



thereby reducing the continuous rotation symmetry to a 
three-fold one. This Z3 freedom arises from the possi- 
bility to translate the configuration of the charge-density 
wave as a whole by a unit cell vector ai or a2. A finite 7 
also shifts the critical interaction strength Vc to a smaller 
value compared to . Moreover, it turns the second or- 
der (quantum) phase transition into a first-order one. A 
crystal-field term cx cos(p(,c) can also strongly affect the 
thermodynamic properties of the model and the value 
of the integer p is important. ^'^ For the planar a;y-model 
supplemented with a crystal-field term cx cos(p(^) it has 
been shown that the ground state always has a broken 
symmetry. However, for p > 4 and at higher temper- 
atures, there is a critical phase characterized by bound 



vortex-antivortex pairs similar to the one found in the 
absence of the crystal-field. They can unbind above the 
Kosterlitz-Thouless-Berezinskii temperature. For p = 3, 
such a critical phase is absent in the planar model. 



IV. Z3-VORTICES 

To study spatially fluctuating solutions we start again 
from the energy functional (18). We find that the gra- 
dient term proportional to 77 appears in the expansion 
because of the Dirac-like single-particle spectrum in mo- 
mentum space at / = 1/3. However, this term disappears 
for slowly varying configurations once a gap is opened in 
the non-interacting system with finite 5t. To keep our 
discussion simple we set 77 = in the following. In this 
case, there are two distinct length scales in the problem. 
One is the coherence length ^ = ^JK/\a\ which describes 
the characteristic length scale over which the amplitude 
of the order parameter changes. The second one is related 
to the anisotropy and naturally appears in the equation 
of motion for Lp: 



1 



siTi{pip). 



(21) 



This is the so-called sine-Gordon equation and in our 
model p — 2>. The characteristic length scale for the 
anisotropy is Ap — ^2k,/ (p7| A|) and it is a sensible quan- 
tity as long as |A| k, const. If the linear extension L of 
the system satisfies L <C Ap the anisotropy term on the 
right hand side can be neglected'^"' and we can consider 
a special class of (singular) solutions to V'^ip — 0: 



•^{x, y) =q Im[log(x + iy)] . 



(22) 



These vortex solutions have an integer nonzero topolog- 
ical charge (vorticity) 



1 

2^ 



Vip ■ ds, 



(23) 



c 



where C is a loop encircling the singularity at the ori- 
gin. The energy of a single vortex configuration grows 
logarithmically with system size L. If ^ ^ L ^ Ap we 
expect that thermally excited vortex-antivortex pairs are 
present. 

On the other hand, on length scales L ^ \p, the right 
hand side of Eq. (21) can no-longer be neglected. Then, 
the simplest nontrivial solution is a domain wall between 
two degenerate ground states. An example describing a 
kink which extends along the x axis is 



(fiy) = \ — atan e 

p p V 



-\/vvl\ 



(24) 



and from the energy functional Eq. (18) it follows that 
there is a finite energy per length associated with the 
domain wall. There exist also single and multi- vortex so- 
lutions of Eq. (21) which are obtained by deforming the 
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vortex solutions of the Laplace equation. (For p = 4 
and g = ±1 a particularly simple explicit expression is 
known.) The single vortex-like solutions have the prop- 
erty that for \R\ <C Xp they reduce to expression (22) 
whereas for \R\ ^ Xp the domain walls between the de- 
generate ground states are resolved, as in Eq. (24). As a 
result, the energy of such a Zp vortex eventually grows 
linearly with system size. 



A. Triangle-rule violation in the vortex center 

Z3 vortices can also be considered in the classical limit 
and two examples are shown in Fig. 4. For a defect 
with (7 = ±1, the complex phase of the order parame- 
ter changes from ip^ to (ps to (pc and back to (pA- This 
situations is sketched in Fig. 4(a). It turns out that for 
such an elementary defect there is necessarily one triangle 
(shaded) where the triangle rule is violated. For an up- 
triangle-rule violation ("up defect"), the phase changes 
clockwise while for a down-triangle-rule violation ( "down 
defect") it changes counterclockwise. Furthermore, as 
explained in Sec. II B, an empty triangle binds a (posi- 
tive) deficit charge of 1 /2 compared to the uniform phase 
while a triangle with two fermions binds —1/2. Thus, 
we can label elementary defects by the pair {Q, S) where 
Q = ±1/2 refers to the bound charge and <5 = A or V 
indicates on which triangle the triangle rule is violated. 
These defects are topologically protected. Note that in 
the classical limit (vanishing hopping), domain walls do 
not cost any energy because the triangle rule is fulfilled. 

It is also possible to construct defects which are com- 
posed of more then one elementary defect. An example is 




FIG. 4. (Color online.) Schematics of two different Z3- 
vortices: A vortex with circulation q — 1 is shown in (a) 
and a vortex with circulation g = — 2 in (b). In (a), three 
domain walls meet at a down triangle and the triangle-rule is 
necessarily violated. As a consequence, a charge 1/2 (—1/2) 
is bound to this topological defect if the center triangle is 
empty (singly occupied). The configuration in (b) shows a 
double vortex where six domain walls meet. In the center 
there is one particle for one down triangle and three up tri- 
angles. This effectively results in two down triangle rule vio- 
lations. The double vortex shown in (b) binds a charge 1 but 
it can be split into two topologically protected single vortices 
with a bound charge of 1/2. 




FIG. 5. (Color online.) Periodically extended charge configu- 
ration of a finite hexagonal system with defects at its corner. 
The considered system contains 1641 sites, 536 particles and 
two well-localized defects (up- and down triangles) each bind- 
ing a charge 1/2. The three different uniform ground states 
(A, B and C) meet at the center of the defects and domain 
walls are indicated by the solid curves. The diameter of the 
circles building the kagome lattice is proportional to the local 
density and the interaction has been set to V = 4t. 



sketched in Fig. 4(b) where the phase changes twice when 
going clockwise around the defect. It can be viewed as a 
composite object of two elementary up defects (1/2, A) 
which in total binds a unit positive charge. Therefore, 
this defect is topologically not protected since it can be 
split into two elementary up defects. Another example 
is the composite object involving (1/2, A) and (1/2, V) 
which can be viewed as a polaron state, see Sec. VC. 



V. NUMERICAL SOLUTIONS WITH DEFECTS 

We now turn to a numerical study of mean-field so- 
lutions with defects at zero temperatures. Thereby, we 
will focus on the properties of the elementary Z3 vor- 
tices as sketched in Fig. 4(a). Because the defects are 
charged, we found it necessary to dope the system in 
order to stabilize mean-field solutions with defects. We 
therefore discuss examples where a single hole has been 
doped into a finite system with periodic boundary con- 
ditions. Self-consistent solutions are found by iterating 
the self-consistency equations. If the interaction is not 
too close to the critical interaction (V > 3t > Vc ^ 2.2t), 
meaning that the defect size is comparable to the lattice 
constant, it is possible to choose the initial charge con- 
figuration such that solutions with two separated defects 
are stabilized. 
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A. Fractionally charged defect 

Let us first look at a configuration wfiere the defects 
form a regular lattice. In the most symmetric case, the up 
and down defects are arranged on interpenetrating trian- 
gular lattices and the defect lattice then has a hexagonal 
symmetry. Such a configuration offers a convenient possi- 
bility to investigate the properties of an isolated defect if 
they are separated far enough from each other. In actual 
calculations we considered finite systems of the form of a 
hexagon with defects located at its corners and employed 
periodic boundary conditions. Figure 5 shows the self- 
consistent charge configuration (periodically extended) of 
such an arrangement. In this example, we have consid- 
ered 1641 sites with 536 particles and the interaction has 
been set to V = At. The diameter of the circles building 
the kagome lattice is proportional to the local density. As 
before, the three inequivalent uniform phases are denoted 
by A, B and C and in Fig. 5 the domain walls between 
them are indicated as solid curves. In the initial state, 
the domain walls are straight lines but after the iteration 
process has converged they reveal a winding character. 

Figure 6 shows the single particle energy spectrum of 
the same system. The defect states have energies which 
lie in the gap of the uniform phase. In total, there are 
six in-gap states which gives three states per defect. Out 
of the three states, one state is lower in energy than the 
other two. For the considered interaction strength, the 
defect size is comparable to the lattice spacing. This 
means that the defect states are basically localized on 
a single triangle and the energy splitting can be de- 
rived from the eigenenergies of a particle hopping on an 
isolated triangle (with amplitude t). Indeed, we have 
checked that in the large V limit the energy splitting be- 
tween the two upper states and the lower one approaches 
Zt — t — {—2t). Note that because the overlap between 
up and down-defects is very small for the considered sys- 
tem, the energy splitting between them is not visible in 
Fig. 6. The step-like features in the spectrum at higher 
energies (near 6t and 7t) are the remains of the subband 
formation due to the enlarged unit cell of the uniform 
charge-density wave. For the uniform solution wc find 
true energy gaps between the steps but in the presence 
of the topological defects, these gaps are filled. The states 
with energies between the steps are localized along the 
domain walls which can be confirmed by studying the 
wave-functions in real space. 

From the analysis of classical charge configurations 
we expect that an elementary defect binds a fractional 
charge ±1/2. We now want to confirm this result for fi- 
nite t by an explicit calculation for the defect lattice of 
Fig. 5. In order to get rid of the short wave-length density 
oscillations, we have considered an averaged charge dis- 
tribution on the hexagonal lattice defined by the center of 
mass points of the triangles forming the kagome lattice, 
as explained in Sec. II B. Figure 7(a) shows the integrated 
density deficit (measured from the uniform particle den- 
sity 1/2) within a circle of radius R around the origin. 
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FIG. 6. Single-particle energy spectrum for tlie defect lattice 
shown in Fig. 5. The inset shows a blow-up of the spec- 
trum near the gap. In total, there are six in-gap states which 
gives three per defect. The spUtting between the upper four 
and lower two in-gap states is proportional to the hopping 
t whereas the splitting among the upper or lower states is 
exponentially small. The quasi continuum of states between 
the step like features at higher energies are the domain-wall 
states. 



The location of the origin has been chosen away from 
a high symmetry point of the defect lattice and is indi- 
cated in panel (b) and (c). As function of R, there are 
clearly visible steps of 1/2 in the integrated deficit den- 
sity which shows that every defect binds 1/2 of charge. 
The steps are rather sharp indicating that the defects are 
rather well localized. This is also seen in panel (b) and 
(c) where the charge deficit and excess measured from 
1/2 in a logarithmic scale is shown. Indeed, most of the 
charge deficit is located very close to the defects but also 
along the domain walls, the density deviates from its uni- 
form value. As a matter of fact, the charge density shows 
oscillatory behavior which is reminiscent of Friedel oscil- 
lations and the density can also exceed its uniform value 
in certain regions, see panel (c). 



B. Irrational versus rational charges 

The value of 1/2 for the charge bound to a topological 
defect depends on a crucial symmetry, namely that the 
up-triangles are equivalent to the down triangles. To see 
this we now consider the effect of a finite 5t in Eq. (2). 
This perturbation effectively increases the hopping on 
the up triangles, = ^ + <5i, thereby breaking the sym- 
metry between the up and the down triangles. We have 
calculated the charge bound to a defect as function of 
dt for different interaction strength. The result is shown 
in Fig. 8. Clearly, the value of the charge varies contin- 
uously with the strength of the symmetry breaking po- 
tential St/t. The effect is larger for smaller values of the 
interaction. This behavior is in agreement with Refs. 13 
and 14 where field-theoretical methods have been used to 
study the effect of a symmetry-breaking potential. Inter- 
estingly, it was found that by introducing a chiral gauge- 
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FIG. 7. Color online, (a) Integrated charge deficit within a 
circle of radius R around the origin. The steps of magnitude 
1/2 indicate that a charge of 1/2 is bound to each defect, (b) 
Logarithmic color plot of the locally averaged charge deficit 
and (c) logarithmic color plot of the locally averaged charge 
excess. In (b) and (c) a cut-off of 10~^ has been used; con- 
sistent with the numerical precision of the solution. The pa- 
rameters of the defect lattice are the same as in Figs. 5 and 
6. 




FIG. 9. (Color online.) The charge distribution of different 
polaron states obtained for V = 3t. In (b) and (d), we have 
marked the misplaced Fermi-rich sites created by separating 
the up and down defects with a cross. 



C. Polaron state 
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FIG. 8. The charge bound to a topological point defect as 
function of the symmetry-breaking field St for interactions 
V = 3t, V = At and V — 5t. The upper set of curves corre- 
sponds to the A-vortices while the lower one corresponds to 
the V-vortices. 



field the energy of a single defect becomes finite turning 
them into well-defined excitations.^^ But if this happens, 
the fractional charge is rerationalized to 1/2.^^ Turning 
this argument around, we may view the dependence of 
the charge on 5t as a manifestation of the fact that in 
our system a chiral gauge field is absent and that a single 
topological point defect costs an energy which depends on 
the system size. We want to stress again that the value of 
1/2 is protected by the symmetry between up and down 
triangles and does not result from a spectral symmetry of 
the single-particle excitations. Such a particle-hole sym- 
metry is absent on the kagome lattice and only emerges 
in the low-energy description, see Sec. VI. 



In the previous section, the geometry and the bound- 
ary conditions have been chosen such that the property 
of an isolated defect can be studied. However, to have a 
configuration which only costs a finite energy in the ther- 
modynamic limit, the vortex has to be healed at some 
point. This is possible when considering pairs of defects 
and naturally shows up when we study the property of a 
single hole doped into a large system. 

In the following, we want to find the ground state 
mean-field solution for a single hole. First, we note 
that the uniform mean-field solution which preserves the 
translational symmetry is always higher in energy than 
several solutions with inhomogeneous charge distribu- 
tions where the inhomogeneity is restricted to a relatively 
small region. This signals the failure of the rigid band 
picture in the interacting system. As a matter of fact, a 
single hole doped into the uniform phase tends to polarize 
its surrounding which leads to an inhomogeneous charge 
distribution around the hole. Following standard nomen- 
clature, we refer to the hole with its polarizing cloud as 
a polaron state. What is special in the present system is 
that the polaron has an internal structure. In fact, it can 
be viewed as a bound state of an up and a down defect. 
The confinement of the two defects results from the en- 
ergy cost of domain walls which are necessarily created 
when trying to separate them in real space. 

Figure 9 shows the self-consistent charge distribution 
of various polaron states for V = 3t. The density dis- 
tribution shown in (a) can be thought as the result of 
removing an electron from a Fermi-rich site. Clearly, 
the polaron wave function is well localized. Panel (c) 
again shows a fairly well-localized polaron state with 
large isotropy. However, in this case the polaron is in 
an excited state. Panel (b) shows the situation where 
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the position of a Fermi-rich site (marked with a cross) 
has changed as compared to (a). As a result, up and 
down defects have been separated and the polaron wave 
function acquires two components. By moving the po- 
sition of other Fermi-rich sites, the two defects can be 
separated even further, as shown in (d). In the sim- 
plest case, this procedure generates a straight "string" 
which connects the two defects. If the charge-density 
wave is in the ground state A at infinity, the string of 
misplaced Fermi-rich sites can be viewed as the phases 
B and C with a minimal extension in the direction per- 
pendicular to the string. (We note that it is also possible 
to stabilize solutions where all the three phases are ex- 
tended but we have found that these configurations have 
higher energy than the straight string for the same sep- 
aration of the two defects.) The energy as a function of 
the number D of misplaced sites [the sites marked with 
a cross in Fig. 9(d)] measures the confinement between 
the two defects. Figure 10 shows the energy of the hole 
as function of D obtained for V — 3t and V — 5t. In 
both cases, the energy grows linearly with D for large 
D. This is in agreement with the expectation that every 
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FIG. 10. The confinement potential between two fractional 
defects in the charge ordered kagome lattice. Shown is the 
energy difi^erence between mean-field states with a single hole 
and the undoped uniform state as function of the number 
D of misplaced sites necessary to separate the two defects, 
see Fig. 9. In the top panel, the nearest neighbor interac- 
tion is y = 3t. The minimal energy occurs at a separation 
D — 3. In the bottom panel, the nearest neighbor interac- 
tion of y = 5t is considered. The minimal energy occurs at 
zero separation of the defects. Results are obtained by solv- 
ing the self-consistency equations on a hexagonal cluster with 
periodic boundary conditions including 1296 sites. 




FIG. 11. (Color online.) The charge distribution of a po- 
laron state in the weakly ordered charge-density wave. The 
interaction has been set to = 2.3t. 



misplaced site costs the same energy because the ring ex- 
change in the hexagons participating in the string is no 
longer effective. The slope is smaller for larger V and 
eventually vanishes for V/t ^ 00. In this limit, the two 
defects are free to separate. For V = 5t, the lowest en- 
ergy configuration is the one with tightly bound defects 
[cf. Fig. 9(a)]. It then costs a finite energy to separate the 
two defects by one unit [cf. Fig. 9(b)] and after that the 
energy increases linearly with distance D. Interestingly, 
the situation looks quite different for weaker interactions. 
Namely, as shown in Fig. 10 for V = 3t, the lowest en- 
ergy configuration corresponds to two defects separated 
by a string of length D = 3. Thus, in this regime, the 
polaron state has a diatomic molecule character where 
the confinement length exceeds the defect size. 

Quantum mechanical processes which go beyond 
the static mean-field description (such as the ring 
exchange^^) will alter the quantitative dependence of the 
confining potential on t/V . However, we expect the qual- 
itative aspects of the static mean-field solutions to be ro- 
bust against a more careful treatment of these processes. 

The above systematic approach only works when the 
defect size is comparable to the lattice spacing which is 
the case for V > 3t. We have also numerically studied 
the polaron wave function for smaller interactions in the 
weakly ordered state close to Vc- A typical converged so- 
lution is shown in Fig. 11 for V — 2.3t. In this regime, the 
polaron extends over several lattice spacings and it was 
not possible to control the location of the elementary up 
and down defects. Moreover, typical solutions are rather 
isotropic. This indicates that the confinement length at 
zero temperature is smaller than the defect size. 

In general, we expect that the polaron is dynamical 
and there is a center of mass motion as well as a rela- 
tive motion of the two defects forming the polaron. For 
a clean system, it is likely that the polaron is not local- 
ized in real space. Rather, one would try to restore the 
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K € Bi (Z = 0, ±) we decompose the crystal momentum 
according to K = Ki+k with k & Bq and define c^i{k) = 
[K) . This allows us to write the local Fermi operators 
as a sum over the patches in the following way 



(25) 



Here, we have separated out the oscillatory factors 
exp[l^; • [r + r^)] and have defined 



FIG. 12. The \/3 x \/3 reconstruction of the charge ordered 
state leads to a partitioning of the Brillouin zone. 



translational symmetry by taking a superposition of lo- 
calized polaron states with the same or nearby energies 
by considering the configuration interaction between the 
different mean-field states. Such an approach has for ex- 
ample been used to study the dispersion of a doped hole 
in the Hubbard model. On the other hand, in the pres- 
ence of imperfections, trapping of the polaron can occur. 

In the strongly correlated regime, the quantum me- 
chanical polaron wave function has a large spatial ex- 
tent because the confining potential is weak. Therefore, 
increasing the doping concentration could lead to new 
quantum phases where the confining is no longer rele- 
vant. For example, one can speculate that a plasma of 
fractionally charged defects is realized once the mean dis- 
tance between polarons fall below the average diameter 
of a single bound pair.'^^ Other interesting phases may 
involve crystalline structures of fractionally charged de- 
fects. 



VI. WEAKLY ORDERED STATE 

To overcome the limitations of the numerical approach 
for the weakly ordered state we now turn to a more an- 
alytical description of topological point defects in the 
regime where both the order parameter and its gra- 
dient are small. Thereby, we are making a connec- 
tion to the solitonic fractionalization mechanism in two 
dimensions. Thus, we assume that we are sufficiently 
close to the phase boundary {V > Vc) such that an ex- 
pansion in the order parameter is justified. In linear or- 
der, only the low energy degrees of freedom in the vicinity 
of the two Dirac points enter. 

For the analytical treatment it is convenient to use a 
notation which is adapted to the ^/i x \/3 reconstruction 
of the unit cell in the ordered state. Hence, we divide 
the first Brillouin zone into three patches Bq and B± 
located around F = — (0, 0) and K± as shown in 
Fig. 12. We will always use the convention that a capital 

or Q is defined in the original Brillouin zone B = ®iBi 
while a small k or q denotes a vector in Bq. Then, for 



feGBo 



(26) 



The states in one patch live on a kagome lattice in real 
space with a threefold enlarged unit cell. The single- 
particle operators in the different patches are related by 
the uniform ordering vector G: 



c^-{k + G) = c^+(fe), 
c^Q{k + G) ^ c^-(fc), 
Ci,+ {k + G) = Ct,o{k). 



(27a) 
(27b) 
(27c) 



Therefore, the Fourier components Pa{Q) with Q close 
to ±G couple the single particle states between different 
patches and most importantly, between the two Dirac 



A. Effective sublattice basis on the kagome lattice 

For the low-energy description it is justified to trun- 
cate the Hilbert space by restricting to the single-particle 
states in the vicinity of the Dirac points at K±, see 
Eq. (15). Therefore, only operators associated with the 
two valleys / = ± are kept. The next step involves a 
fe-independent transformation from the site to the "sub- 
lattice" basis: 



ci,±(fc) 



C2,±(fc) 



C3,±(fc) 



V3 



a± 



(fc) -fe±*'^/^6±(fc)l , (28a) 



e±"/3a± (k) + e^'^'^Hi (fc)! , (28b) 



[a±{k) + b±{k)]. 



(28c) 



Above, we have suppressed the contribution of operators 
acting on states of the flat band at higher energy. The a 
and the b operators are chosen in analogy to graphene in 
which case they would act on states living either on the 
A or the B sublattice. It turns out that on the kagome 
lattice, the up and the down triangles play the role of 
the A and B sites. This becomes clear when inverting 
the relation (28) for 0^(0) and bj. (0) and expanding in 
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terms of real space operators:^° 
4(0) 



aUo) 



r 
r 

[wc\{r)+uj^clir) + cl{r) 



/3N 



t fri\ — 



61(0) 



7^ E [4(0 + 4(0+4(0 



(29a) 
(29b) 

■r 

7 

(29c) 
(29d) 



Here, we have introduced a; = exp(27ri/3) and have 
set = for clarity. In the state created by a\^{0) 

[aL{0)], the phase on every up triangle increases by 27r/3 
along each bond in the [anti-] clockwise direction while 
the phase remains constant on the down triangles. On 
the other hand, in the state created by b^_{0) [^^^(O)], 
the phase on every down triangle increases by 27r/3 
along each bond in the [anti-] clockwise direction while 
the phase remains constant on the up triangles. 



B. Projected mean-field Hamiltonian 

The operators introduced in the previous section al- 
lows one to obtain an effective low-energy Hamiltonian. 
The calculation is straight forward but lengthy. In the 
following we will only present the final results. 



1. Kinetic energy 

Linearizing in k around K± and applying the transfor- 
mation (28) brings the low energy tight-binding Hamil- 
tonian into the canonical Dirac form 



Hq = YY1 [^Fl{ka: + ilky)a\{k)bi{k) + h.( 



(30) 



l=± keBo 



with the Fermi velocity vp = {h = 1). We have 

shifted the zero of energy to the Dirac points. 



2. Bond order term 

The term Eq. (2) which breaks the symmetry between 
the up and down triangles enters the low energy descrip- 
tion as a staggered potential fig — —3/2dt for the effective 
sublattice states: 

Hbow - -M. E E HC^)"^^) - b]{k)bi{k)\ . (31) 



l=± keB, 



Here, terms of order 0{k'^) have been neglected. The 
above form is in full analogy with a staggered potential 
on the honeycomb lattice. Furthermore, the effect of a 
finite fig is plausible when considering the representation 
Eq. (29). 



3. Mean-field interaction 

In order to project this operator onto the low-energy 
degrees we first write the density operators n^{—Q) in 
terms of the patch operators c^iik). Then, using the 
transformation (28) and keeping only operators which 
act on the low-energy degrees alone, the mean-field in- 
teraction assumes the following form 



Hy ^Y<^\k - q)V{q)<P{k) + const. 



(32) 



k,q 



Here, the summation is over the reduced Brillouin zone, 
k and q G Bq. Furthermore, we have introduced the 
four-component spinor 



$(fc) 



a+{k) 
a_(fc) 
\b-ik)J 



-y<f>(i?)e- 



ik-R 



(33) 



R 



Note that ^{R) is a coarse grained operator defined on 
a triangular lattice with a three-fold bigger unit cell. We 
assume that both |A| and the gradient V^p are small and 
we keep only terms entering linear in these quantities. 
For the matrix V{q) we then find the following simple 
expression 



V{q) 



A(q)l 
A*(g)l ^ 



(34) 



Here, 1 is the 2x2 identity matrix and 
2V 



A(g) 



3iV 



[pi(G + q)+p2{G + q)+ P3{G + q) 



(35) 



is assumed to be peaked around q « V(/3. 



4. Continuum limit 

The continuum limit is defined by a — >■ while keeping 
the Fermi velocity vp and the coupling constant V ~ 
Vvuc constant. It then follows that we scale the fermion 
fields according to 



WucE ^ I (fR..., $^* = $/, 

R •' 



(36) 



where the unit cell volume is Wuc = 3\/3a^/2. The lin- 
earized mean-field theory takes the following continuum 
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form 



Here, n is a unit vector defined as 



^\R)tC^{R) + —\/^{R)\ 



(37) 



for the four-component wave function ^'(-R) given in 
Eq. (36). In the notation similar to Ref. 12 the kernel 
K, in Eq. (37) is written as 

lC = vpa- {~iV) + + P [Ai (fi) - ij5^2{R)] ■ (38) 

We used the 4x4 Dirac matrices 

(39) 

where i = x,y and the 2x2 matrices and R are defined 
as 



1 



-t\ , p /I 

cr,, = and H — Or = 



10/' \i I " \0 -1, 



Ai and A2 are real and imaginary part of A. 



(40) 




(43) 



The total charge bound to the vortex can be obtained by 
integrating the density over space 



Q = y dxdyp{x,y). 



(44) 



From Eq. (42) it follows that the integral Eq. (44) mea- 
sures the area (in units of Att) which is covered on the 
sphere by the unit vector n in the mapping {x,y) i— > 
n{x,y). The result is therefore 



sign(^s) - 



(45) 



In the limit fig — > 0^ and for g = ±1 we recover the value 
Q = ±1/2. On the other hand, the dependence of the 
charge on fig is similar to the result reported in Fig. 8. 



C. Vortex solution and mid gap states 

We are now in a position to study the effect of a vortex 
in A(ii) on the fermionic spectrum. Thereby, we apply 
the results previously found for graphene.^^'^'^'^'*''*^ We 
assume a symmetric vortex configuration with vorticity 
q which in polar coordinates is written as 

A(r, 0) = Ai+ iA2 = Ao(r)e"?^+". (41) 

The amplitude Ao(r) vanishes for r — > and assumes a 
constant value Ao(oo) far away from the origin. It has 
been shown that such a configuration leads to a single 
mid-gap state at an energy E — fig or E = —fig, depend- 
ing on the sign of the vorticity 5.^^ These solutions merge 
into a single zero-energy mode in the limit /is — >■ 0. From 
the emergent spectral symmetry of the Dirac equation 
and the completeness of states in the absence and pres- 
ence of a vortex it follows that both the valence and the 
conduction band transfer half a state to the zero-energy 
mode. As a result, a charge ±1/2 is bound to the vor- 
tex depending on if the zero-energy mode is occupied or 
not. These simple considerations break down for a finite 
Us because the emergent particle-hole symmetry of the 
single particle spectrum is violated. The calculation of 
the bound charge in this situation is more involved. We 
follow here the argumentation of Ref. 13 and 14 which 
makes use of the fact that a inhomogeneous static con- 
figuration of the three real fields /is, Ai and A2 induce 
a fermionic charge density. This charge density is ob- 
tained from a perturbative treatment around the uniform 
solution:^*' 

p{x, y) = -^n • {dxfi X dyu). (42) 



VII. CONCLUSIONS 

We have studied topological point defects in the 
charge-density wave realized in a model of interacting 
spin-polarized fermions on the kagome lattice at filling 
fraction 1/3. We have found that elementary point de- 
fects carry a charge ±1/2 as long as the symmetry be- 
tween the up and the down triangles of the kagome lat- 
tice is preserved. If this symmetry is violated, the bound 
charge varies continuously with the symmetry breaking 
term. Moreover, we have argued that in the classical 
limit the point defect corresponds to a local violation 
of the triangle rule and is therefore related to the fact 
that the interaction is frustrated. On the other hand, in 
the weakly ordered state, we made a connection to the 
solitonic fractionalization mechanism based on the Dirac 
equation with a vortex in the background field. The con- 
sidered system therefore offers a unique possibility to re- 
alize these two different routes to fractionalization in the 
same model. 

Using unrestricted mean-field calculations we have 
studied the ground state of a single hole doped into the 
charge-density wave. We have found that the polaron 
state can be viewed as a bound state of two defects both 
carrying a charge 1/2. We have also calculated the con- 
fining potential between these two defects and have found 
that in the intermediate interaction regime it is mini- 
mized for a finite separation. 

We note here that the charge ordered state considered 
in this article bears some similarity with the "trimerized" 
phase considered previously in that it shares the same 
enlarged unit cell and also has three different ground 
states. ^^'^'^ However, the charge density- wave order seems 
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to be more easily realized in an interacting system. Our 
basic conclusions remain valid when the spin degree of 
freedom is taken into account as long as the interaction 
still favors the charge-density wave. Nevertheless, in the 
spinfull case it is not the charge which fractionalizes but 
the defects carry either a spin 1/2 and no charge or a 
charge ±1 and no spin. 
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